if !d.name eq 'PS' then begin
    device,xsize=12,ysize=12,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end
@eu_setup
!p.multi=[0,1,4] & !p.charsize=2.
nt=(size(bx))[3]
print,nt
nl=64
year = 365.*24.*60.*60.
time = (nxaver*dt*indgen(nt))/year
good=where(time gt 40)
ngood=(size(good))[1]
bmax=max(abs(bx[*,*,good]))
bmin=-bmax
lev=grange(bmin,bmax,nl)
if (itac eq 0) then begin
 rlev=33
endif else begin
 rlev=62
endelse
bphi=reform(bx[*,rlev,good])
bphi_r=reform(bx[25,*,good])
bth=reform(by[*,rlev,good])
br=reform(bz[*,rlev,good])

lev = 2.*max(abs(bphi))*indgen(nl)/float(nl-1) $
        -max(abs(bphi))

cpos=[0.13,0.72,0.81,0.95]
bpos=[0.95,0.72,0.96,0.95]

bphi_or=fltarr(l,m,ngood)

; organizing the magnetic field
for k=0, l-1 do begin
  for j=0,m-1 do begin
    bphi_or[k,j,good] = bx[j,k,good]
  endfor
endfor

contour,rotate(bphi,3),time[good],(th),/fi,nl=64,ytitle='!8B!D!7u!N!6',$
;contour,transpose(reform(bphi_or[rlev,*,good])),/fi,nl=64,ytitle='!8B!D!7u!N!6',$
        ystyle=1,lev=lev,pos=cpos,xstyle=1
colorbar,range=[min(lev),max(lev)],/vertical,xaxis=2,$
         ytickformat='(F6.2)',yticks=4,pos=bpos,charsize=1.5,$
         ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)]


cpos=[0.13,0.42,0.81,0.65]
bpos=[0.95,0.42,0.96,0.65]
lev = 2.*max(abs(bth))*indgen(nl)/float(nl-1) $
        -max(abs(bth))
;contour,rotate(bth,3),time[good],th,/fi,nl=64,ytitle='!8B!D!7u!N!6',$
;        ystyle=1,lev=lev,xstyle=1,pos=cpos
lev=grange(-1.0,1.0,nl)
contour,smooth(rotate(bphi_r,4),2),time[good],r,/fi,nl=64,ytitle='!8B!D!7u!N!6',$
        ystyle=1,lev=lev,xstyle=1,pos=cpos
colorbar,range=[min(lev),max(lev)],/vertical,xaxis=2,$
         ytickformat='(F6.2)',yticks=4,pos=bpos,charsize=1.5,$
         ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)]


cpos=[0.13,0.12,0.81,0.35]
bpos=[0.95,0.12,0.96,0.35]
lev = 2.*max(abs(br[5:59,*]))*indgen(nl)/float(nl-1) $
        -max(abs(br[5:59,*]))
contour,rotate(br,3),time[good],th,/fi,nl=64,ytitle='!8B!Dr!N!6',$
        xtitle='time [years]!6',ystyle=1,lev=lev,xstyle=1,pos=cpos
colorbar,range=[min(lev),max(lev)],/vertical,xaxis=2,$
         ytickformat='(e11.2)',yticks=4,pos=bpos,charsize=1.5,$
         ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)]
time=time[good]
save,bphi,bth,br,th,time,filename='bfield.sav'

!p.multi=0
end
